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Construction of power flow feasibility sets 

Krishnamurthy Dvijotham, Konstantin Turitsyn 


Abstract —We develop a new approach for construction of 
convex analytically simple regions where the AC power flow equa¬ 
tions are guaranteed to have a feasible solutions. Construction 
of these regions is based on efficient semidelinite programming 
techniques accelerated via sparsity exploiting algorithms. Result¬ 
ing regions have a simple geometric shape in the space of power 
injections (polytope or ellipsoid) and can be efficiently used for 
assessment of system security in the presence of uncertainty. 
Efficiency and tightness of the approach is validated on a number 
of test networks. 

I. Introduction 

Future power systems relying on large amounts of clean and 
renewable generation will operate at highly increased levels 
of uncertainty. Lack of controllability and predictability of 
renewable generation output will require substantial revision of 
modern planning and operational practices. On the operational 
level, unexpected variations of wind and solar power can po¬ 
tentially compromise the otherwise secure systems. From the 
long-term perspective, heavy reliance on weather dependent 
distributed generation sources increases the spatial variability 
of load levels, and compromises the validity of approaches 
based on scenario analysis. 

In recent years there has been an explosion in the num¬ 
ber of academic works that extend the existing procedures 
by incorporating uncertainty. Most common approaches in¬ 
clude stochastic programming, robust and chance-constrained 
counter-parts of traditional optimal power flow problem Q- 
Q. Despite the progress achieved by the academic com¬ 
munity, the proposed algorithms are usually computationally 
prohibitive and may not be adopted by industry in the near 
future. Moreover, most of the existing approaches rely on 
linearized versions of power flow equations and may not be 
suitable for operational security assessment purposes where 
nonlinear effects dominate. 

Our work attempts to address the need for computationally 
tractable tools to assess the effect of renewable uncertainty on 
power system security. The key contribution of our work is a 
technique for characterization of a maximal set of uncertain 
power injections that can be tolerated by the power system 
while maintaining feasibility. We consider two types of un¬ 
certainty sets: first is a polytope, where the power injection 
on a number of buses can vary independently. The second is 
an ellipsoid which represents chance constraints in situations 
where variability of individual resources is correlated and has 
a Gaussian probability distribution. For both of the sets our 
construction guarantees that the AC power flow equations will 
have feasible solution for any injections in the uncertainty set. 
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Moreover, by construction the set is a maximal for a given 
shape. In other words any uniform rescaling of the set will 
result in violation of one of the constraints or disappearance 
of the solution. 

Characterization of the feasibility region has been addressed 
by multiple authors in the recent years. Several studies in the 
last decades have looked at the question of convexity of the 
feasibility set jbl- Non-convexity of the feasibility region 
for general networks has been explicitly demonstrated in Q, 
0 - In the last decade, however, much progress has been 
made in understanding sufficient conditions for convexity of 
the feasibility region 0-0 relevant for global optimality 
of optimal power flow relaxations. Under the assumption that 
the network is lossless and all buses are PV buses, a number 
of papers have studied necessary and sufficient conditions that 
guarantee the existence of power flow solutions 0- In recent 
work 0^ the authors propose a general framework to con¬ 
struct Linear Matrix Inequality (LMI)-based inner approxima¬ 
tions to the feasible set of a class of quadratically constrained 
quadratic programs (QCQPs). However, the framework does 
not deal with arbitrary quadratic equality constraints like the 
AC power flow equations for a meshed network. 

Computationally tractable characterization of the power 
flow feasibility region in terms of ellipsoids in injection space 
was proposed in 0. This is similar to the regions developed 
in our paper, but relies on the linearized DC power flow 
approximation. A dual characterization of the feasibility set in 
terms certificates of insolvability has been developed recently 
in 0 - More recently, a do-not-exceed limit strategy was 
proposed for optimal dispatch instructions for intermittent 
generation p^ , pT] , that relies on similar mathematical 
constructions as developed in this paper but limits the analysis 
only to DC power flow models. 

A concepmally close paper that has largely inspired our 
effort is the recently published study 0 presenting the 
construction of certificates for existence of AC power flow 
solutions. This work relies on well-known Banach fixed 
point theorem to construct ellipsoidal and polytopic regions 
where the power flow solutions are guaranteed to exist. This 
has been extended in 0 to reduce the conservativeness 
of the constructed certificates. This manuscript extends and 
generalizes the previous approaches, and proposes a unified 
methodology for constructing the certificates for existence 
of feasible solution to power flow equations. In comparison 
to previous studies it allows for natural incorporation of all 
the important constraints on voltage levels and power flows, 
and at the same time provides a computationally tractable 
methodology based on semideflnite Programming (SDP) that 
greatly improves the size of the certified region in comparison 
to previous studies. 

Presentation of our results is organized as follows. In section 
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|I^we introduce the notations and key modeling assumptions 
employed in the paper. Section III introduces the mathematical 
background behind our approach. Applications of the approach 
in practice are described in IV In section |V] we validate the 
approach via analysis of several standard IEEE cases. We 
conclude in section [Vl| by assessing the results and discussing 
possible extensions of the approach. 


II. Modeling 

In this section we define the mathematical notations used 
throughout the paper, and discuss the modeling assumptions 
behind the construction of certified feasible regions. 


A. Notations 

We use M to denote the set of real numbers, C the set of 
complex numbers. M", C" denote the corresponding Euclidean 
space in n dimensions. Given a set C C M", Int (C) denotes 
the interior of the set. Given a complex number x G C, 
Re (x) denotes its real part and Im (x) its imaginary part. 
X* denotes its complex conjugate. |a;| is its absolute value and 
Zx G [— TT, tt) is its phase. 1 denotes the vector with all entries 
equal to 1. The Jacobian of a function F : R" i—> K™, denoted 
as Jf, is an n X TO matrix whose i-th row is the gradient of 
Fi. Given a finite set S, 151 denotes the number of elements 
in the set. 


B. AC power flow model 

We represent the transmission network as a graph (V,5) 
where V is the set of nodes and 8 is the set of edges. In 
power systems terminology, the nodes represent the buses and 
the edges correspond to power lines. Buses are denoted by 
indices i = 0,1,... ,n and lines by ordered pairs of nodes 
(i, j). We pick an arbitrary orientation for each edge, so that 
for an edge between i and j, only one of (i, j) and (j, i) is in 
5. 

The transmission network is characterized by its complex 
admittance matrix Y G C"^". Y is symmetric but not 
necessarily Hermitian. Define G = Re (F), i? = Im (F). 

Let Vi be the voltage phasor, pi and qi denote active and 
reactive injection at the bus i respectively. V is the vector 
of voltage phasors at all buses. Three types of buses are 
considered in this work; 

• PV buses where active power injection and voltage mag¬ 
nitude are fixed, while voltage phase and reactive power 
are variables. The set of PV buses is denoted by pv. The 
voltage magnitude set-point at bus i S pv is denoted by 

Vi. 

• PQ buses where active and reactive power injections are 
fixed, while voltage phase and magnitude are variables. 
The set of PQ buses is denoted by pq. 

• Slack bus , a reference bus at which the voltage magnitude 
and phase are fixed, and the active and reactive power 
injections are free variables. We choose bus 0 as the slack 
bus as a convention. 

We denote the union of PV and PQ buses as nsb = pv U pq. 
Let k = |nsb| -f |pq|. This is the total number of variables to 
be solved for in the power flow equations. 


Definition 1 (Valid Voltage Phasor Vector). A vector V £ 
C”+\ indexed by i — 0,..., n, is said to be a valid 
voltage phasor vector if |Fi| = Vi for each i £ pv and 
|K)| = vq,ZVq = 0. Throughout this paper, we will work 
only with valid voltage phasors. We denote the constraints on 
valid voltage phasors as (F) = 0. 

Definition 2 (Injection Vector). A vector s £ is said to 
be an injection vector. The first n coordinates coiTespond to 
active power injections at buses 1 through n, and the last 
|pq| components correspond to the reactive injections at the 
PQ buses. 

Definition 3 (Power Plow Operator). Let F be a valid voltage 
phasor with Vi = exp{pi+j0i). Define the power flow 
operator F as 


n 

[F (F)]j = ^ Bij exp {pi + pj) sin {Oi - 9j) 

j=o 

n 

+ X! {9^-9j),i = l,...,n (la) 

j=o 

n 

[F(F)W. = E Gij exp {pi -f pj) sin {9i - 9j) 
j=o 

71 

- ^ exp (pj -f pj) cos {9i-9j),i = l,..., |pq| (lb) 

j=o 


Definition 4 (Power Plow Equations). Por any valid voltage 
phasor vector F and any injection vector s, the power flow 
equations are give by F (F) = s. The variables solved for are 

/ ^nsb j denote by Jp (V) the Jacobian of F with respect 

VPpq/ 

to 

VPpq/ 

We now quote a result from our related paper pO) that 
expresses the Jacobian matrix as a quadratic function of the 
vector of voltage phasors F. We will use this result in the 
certification procedure described in section |III-A 

lemma 1. The power flow Jacobian Jp (F) can be written as 
a quadratic matrix function of the voltage phasors: 


^A,|F,p+ ^ ryRe(FF/) + ^,,Im(F,F/) (2) 

iGpq 


where A^, A>ij are kx k matrices that are functions of the 
network admittance matrix Y. 

1) Operational Constraints on Voltages: Apart from the 
power flow equations that describe the conservation laws and 
behavior of generators and loads we also consider the opera¬ 
tional constraints (line limits, bounds on voltage magnitudes 
etc.). In particular, we require that the solutions of power flow 
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equations satisfy the following constraints: 



(3a) 

\Vk, - < /fc,/c = 1,...,TO 

(3b) 

< lm{V^{YV)*) <q„iG pv 

(3c) 

q,<yn{VoiYV)l) <q. 

(3d) 

P^<Re{VoiYV%) <Po 

(3e) 


The hrst two constraints can be interpreted as stan¬ 

dard operational constraints imposed on voltage magnitudes 
and scaled current flows in a power system. Additionally, we 
also require the constraints arising from the specification of the 
power flow problem: upper and lower bounds on the reactive 
power injections at the PV buses (generator reactive limits): 
( |3^ , upper and Lower bounds on the reactive power injections 
at the slack bus ( [3dl i and upper and lower bounds on the active 
power injection at the slack bus ( |^ . These constraints are 
collectively denoted as H°p (T^) > 0 

Although only these constraints are enforced in our simula¬ 
tions, the technique can be naturally extended to incorporate 
any constraints defined via quadratic inequalities in V. 

III. Algorithmic construction of feasible region 

Our goal in this paper is to And regions in the space of 
injections such that every point in the region has a power 
flow solution satisfying the operational constraints ([^. We first 
define this precisely: 

Definition 5 (Domain of Strict Feasibility (DSF)). s is said 
to be strictly feasible (SF) if there is a solution of the power 
flow equations F {V) = s such that F[°p (V) > 0. We say that 
5 C is domain of strict feasibility (DSF) if every s G 5 is 
strictly feasible. 

Remark 1. We require strict feasibility to deal with mathe¬ 
matical issues that arise in our certification procedure. From 
a practical perspective, this does not matter much, since this 
simply corresponds to make the bounds (line limits, voltage 
magnitude limits etc.) tighter by a very small amount. 

In this paper, we will assume that S is an ellipsoidal or box 
region centered at around a nominal strictly feasible s°. These 
may be the 0 injection vector or a nominal operating point 
known to the system operator to be strictly feasible. We will 
attempt to And the largest ball around that is a DSF. The 
construction is algorithmic and based on the following idea: 
Let s G 5 be arbitrary. We consider the straight line path 
from s° to s. Let be a solution of F (V) = s such that 
fjop (y°) > 0. If Jf is non-singular, the power flow 
operator is locally invertible (by the inverse function theorem) 
in a small ball around s°. Thus, we can take a small (but 
finite) step from towards s, on the line segment connecting 
them. Let us say the new point is s^. Suppose that is 
strictly feasible as well, and = F (l/^) ,F[°p > 0. 

If Jp is non-singular, we can again move a small step 
towards s to get s^. This procedure can continue as long as 
we never hit a singularity of the Jacobian Jp {V) or a point 
such that F[°P (V) f 0. The following sections illustrate how 
we can check this condition using semidefinite programming. 



Fig. 2: Tightening the operational constraints to ensure that 
H°P (V; 7 ) > 0 ^ det [Jp (V)) f 0. 


The procedure is essentially like the continuation power flow 
algorithm m- However, our goal is not to find a power flow 
solution (like in the continuation power flow), but to certify 
that a strictly feasible power flow solution exists for every 
s G S. This is illustrated pictorially in figure 

A. Mathematical Characterization of feasible region 

We now formalize the arguments outlined in the previous 
section. We start by making an assumption on S: 

Assumption 1. S is convex and 3s° G S that is strictly feasible. 

Remark 2. In most cases, s° = 0 satisfies these conditions so 
this assumption is not problematic. More generally, we will 
attempt to certify that all injection vectors in a ball around 
a known strictly feasible nominal injection vector are strictly 
feasible. 

Theorem III.l. Suppose the following implications hold: 

H°P {V) > 0, (y) = 0 ^ det {Jp {V)) f 0 (4a) 

{V) > 0, {V) =0,F{V)gS H°P {V) > 0 

(4b) 

Then, S is a DSF 

Proof: See Appendix. ■ 


B. Checking the Conditions Using semidefinite programming 

Written in terms of the vector of voltage phasors, these 
constraints can be expressed as polynomial equalities or in¬ 
equalities in V. These constraints are generally non-convex 
constraints and may not be amenable to tractable computation. 
Thus, we attack the problem using the moment relaxation 
approach p^ , leading to a semidefinite programming based 
sufficient condition for the implication We do this in two 
steps: We first find a tightening of the operational constraints 
such that every point in the feasible set defined by the 
operational constraint also has a non-zero power flow Jacobian 
- this certifies ( |4al i. Then, we show that if s G 5 is feasible 
injection vector, then it must be strictly feasible - this certifies 

(i§. 
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Voltage Space Injection Space 



Fig. 1; Pictorial depiction of our Framework; The red curve is a closed surface in the voltage space defined by det {Jp {V)) = 0 
and the inner brown curve represents the boundary of the operational constraints on voltages. Inside this region, voltages satisfy 
H°P (]/) > 0. The blue ellipse denotes the region S in injection space. As we move from to s in injection space, the solution 
moves from to V in voltage space along the curve. Our certification procedure consists of ensuring that for any s G S, the 
curve in the voltage space never intersects the red curve or the boundary of the operational constraints H°p {V) = 0. 


1) Non-Singularity of the Jacobian: Our first step is to 
ensure that every valid voltage phasor vector V satisfying 
the operational constraints H°p (V) > 0 (and the validity 
constraint [V) = 0) also satisfies det {Jp (1^)) ^ 0. 
However, this may not be true. Thus, we also allow for a 
tightening of the operational constraints such that this becomes 
true (see figure |^. Specifically, we consider tightening the line 
limits: We replace | — Vjj < fij by |L - Vj I < min ( 7 , fi^) 

for some 7 > 0. Denote the modified operational constraints 
as H°P (Id; 7 ) > 0. We then solve the following problem: 


Maximize 7 such that 

H°P {V; 7 ) > 0, {V) = 0 det {Jp {V)) 7 ^ 0 (5) 

The implication in is equivalent to infeasibility of the 
following system: 


Find {z,V) such that ( 6 a) 

Jp {V)z = 0 ( 6 b) 

z'^z = 1 ( 6 c) 

iF°P (y; 7 ) > 0 , (Id) = 0 ( 6 d) 


If this system is infeasible, we know that for each V such 
that H°P (y-jy) > = 0 , there is no z such that 

z^z = 1 and Jp {V) z = 0, which means that Jp (Id) is 
non-singular. This feasibility problem involves non-convex 
constraints and is difficult to solve in general. However, 
all the constraints are defined by polynomial equations and 
inequalities in (V) z). We perform a moment relaxation of 
this problem p^ , which replaces the non-convex constraints 
by weaker convex constraints. If there is no (z, Id) that 
satisfy the weaker constraints, then we know that (|^ is also 
infeasible. Checking the infeasibility of the moment relaxation 
to solving a semidefinite program, for which there are well- 
known efficient algorithms and software. This conversion is 
discussed in Appendix section VII-B We note the sizes of the 
resulting SDPs for some IEEE benchmarks in section |V] 


2 ) Strict Feasibility of Constraints: Eor each constraints i 
in iT°P, we solve the following feasibility problem; 

Eind Id such that (7a) 

H°P(1/;7)=0 (7b) 

H°P (E; 7 ) > 0, (E) = 0, F’ (E) = s, s G 5 (7c) 

If the above problem is infeasible, then we have a certificate of 
( |4b| ). The constraints are quadratic functions of E. Again, this 
is a non-convex optimization problem, but it can be relaxed 
using the moment relaxation approach, described in section 
|VII-B| The overall certification procedure is summarized in 
Algorithm [T] 

3) Tightness of the Certificate: In practical cases (section 

0 , we find that the value of 7 obtained by solving © often 
exceeds the practical bounds by a significant margin. Thus, 
for all practical purposes, we can assume that the tightening 
in step 1 of the certification procedure does not significantly 
impact the conservatism in the certificate. Now let us consider 
the second part, which involves checking infeasibility of Q 
for each i. Suppose that this fails for some i. Do we know 
that S is not a domain of strict feasibility? In other words, can 
one be assured that there exists s G S such that no solution of 
F (E) = s satisfies H°p (E; 7 ) 0? In general, the answer is 

no (if it were not, then would be solving an NP-hard problem 
exactly). However, when the infeasibility test fails, the solver 
does output (E, s) which are feasible for a relaxed version of 
the constraints in Q (see section [VII-B| i. One can then check 
if (E, s) satisfy the constraints in the problem (0. If this is 
true, one is certain that the set S is not a domain of strict 
feasibility. In the numerical section |V] we find that this is 
indeed the case for the test cases we experiment with. 


Algorithm 1 Certifcation Procedure 

Eind maximum 7 such that is infeasible. 

Check that the relaxation of Q is infeasible for each i. 
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IV. Practical applications 

Despite the promising opportunities for acceleration of the 
SDP solvers described in the previous sections, the resulting 
problem is still going to be too computationally intensive 
for real-time applications. Hence, we propose to use the 
algorithm in two stage procedure. On the first, offline stage 
one or multiple certificates are constructed for different base 
operating points V^, and possibly different topologies of the 
system. Each of the certificates establishes feasibility of the 
whole region of injection space, and can be reused multiple 
times in online applications. 

During the online stage the simple algebraic form of the 
region definition is used for fast assessment of feasibility of 
possible operating points and for fast convex optimization of 
control actions. We envision several applications where the 
certificates could be used in real-time operation: 

Security assessment. Many of the real-life contingencies can 
be represented as rapid changes in power injections. These 
mainly result from triggering of protection systems on loads 
and distributed generation, although in the future events like 
unexpected wind gusts and clouds covering major solar plants 
may have similar effect on the system. Certificates established 
in this work can be used to screen the safe scenarios from 
the contingency list, and identify critical contingencies. More 
generally they naturally establish the levels of uncertainty in 
power injection vector that can be tolerated by the power 
system. Whenever the prior distribution of the renewable 
fluctuations is known, the certificates can be also used to 
verify the low probabilities of system failure. This can be 
accomplished by inscription of the injection probability level 
set in the certified feasibility region. 

Establishment of do-not-exceed limits. As proposed in lEZl 
the do-not-exceed (DNE) limit strategy is viable option for 
limiting the effect of uncontrollable renewables on the security 
of the system. The certificates can be used for regular updates 
of DNE limits. Although this strategy can be more conser¬ 
vative in comparison to optimization approaches developed 
in 1161, its advantages include reliance on AC power flow 


equations and effectively zero computational cost. The latter 
advantages opens the possibility of using the approaches in 
micro-grid setting with limited computational resources for 
dispatch decision making. 


V. Simulations 

As noted previously, the first part of the certification ^ only 
depends on the network structure and not on the particular 
operating point. Thus, as long as the network topology and 
parameters are fixed, this can be computed once (offline). In 
table |I] we have the maximum value of 7 (a uniform upper 
bound on | V) — Vj | for all transmission lines) for three different 
test networks. The values obtained show that the bounds are 
non-conservative and looser than the actual flow constraints 
imposed in practice. In the subsequent sections, we apply our 
methodology to three networks included with the Matpower 
package | [23| . We construct polytopic and ellipsoidal regions 
centered around a nominal solution, and check whether the 
inscribed regions are tight (in the sense that they cannot be ex¬ 
panded without violating some of the operational constraints). 


A. 3-bus system 

We start our discussion with a three-bus system model 
depicted in figure]^ The example is taken from p4| and was 
used to prove that the SDP relaxation of the Optimal Power 
Plow (OPE) p5| problem may not be exact . 

All the buses are PV buses, so the only controllable in¬ 
jections are the active injections at buses 2 and 3 (bus 1 is 
the slack bus). The true feasible set (obtained via brute-force 
discretization) is plotted (in the P 2 -P 3 space, in p.u system with 
a lOOMVA base) in blue in figures |3b|3c Superposed on it, 
we plot various subsets S that we certify using the procedure 
from algorithm [T] In figure we look at S defined in terms 
of individual bounds on p 2 , ps and plot the largest rectangles 
with a given center and aspect ratio that can be embedded in 
the feasible set. The results show that the embedding is tight, 
ie, one of the corners of the rectangle is typically close to the 
boundary of the feasible set so that the rectangle cannot be 
expanded while still being a DSP. 

In figure we plot an ellipse inscribed in the feasible 
set. The ellipse covers a much larger area relative to the 
rectangles, since the axes of the ellipse are aligned roughly 
with the directions in which the feasible set extends. Again, 
we see that the ellipsoid is non-conservative in the sense that 
on increasing its radius (while maintaining the aspect ratio), 
it would no longer be contained within the feasible set. 

The ellipsoid embedding also shows how this may be 
advantageous compared to traditional SDP relaxations of the 
OPE problem | [25| p6) . The standard SDP relaxation of the 
OPE problem would work with the convex hull of the blue 
feasible set, which is denoted by the black curve in figure 
In the top left corner of the blue region, the convex hull 
contains points not in the original feasible set. If p 2 is fixed 
to —.6, and the OPE objective is to minimize ps, it is easy 
to see that the SDP relaxation of the OPE problem will find 
a solution that is infeasible for the original problem. On the 
other hand, with our approach, by first inscribing an ellipse 
and minimizing p^ subject to the constraint that {p 2 ,P 3 ) lie 
in the ellipse, we will obtain a point on the boundary of the 
ellipse that is feasible and near optimal (since the boundary 
of the ellipse in that direction is close to the boundary of the 
feasible set). 


B. 6 -bus system 

Next, we consider the 6 bus system shown in figure 
Buses 1-3 are PV buses and buses 4-6 are PQ buses. We 
choose a uniform bound on the voltage differences between 
neighbors \Vi — Vj\ < .4. is indeed infeasible with 
this choice of 7, so that every voltage vector satisfying 
these constraint has a non-singular Jacobian. The injection 
parameters are the active injections at buses 2 and 3, and the 
active and reactive injections at buses 4,5 and 6, leading to 
a 8-dimensional parameter space. We attempt to embed an a 
hypercube (5 = {s : ||s||^ < (5}) inscribed in the feasible set. 
This is the analog of the rectangular region in the previous 
case. We are not limited to this set, but we chose this for 
simplicity. 
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(a) 3-bus network 



(b) Embedding rectangles in feasible set: Blue 
region is the true feasible region and each black 
box is an embedded rectangle in the feasible 
region. 

15 


05 


•05 


•15 


•25 i * i i. .fc t 

•1.S -1 -0.5 0 05 1 15 

(c) Embedding ellipses in feasible set: Blue 
region is feasible region and black ellipse is 
embedded in the feasible region. The hlack curve 
denotes the convex hull of the feasible set. 

Fig. 3: 3 Bus Network 



We compute the maximum S such that S remains contained 
in the feasible region- it turns out to be .728 p.u. This means 
that if all active (at the non-slack buses) and reactive injections 
(at the PQ buses) are bounded by .729 (in absolute value), 
we are guaranteed that there is a PF solution satisfying the 
operational constraints. On increasing 5 further (to .73), we 
find that the upper bound on |Vi — V 5 I is violated (that is 0 
fails to remain infeasible for this constraint). 


System 

Case 3 

Case 6 

Case 14 

Maximum 7 

1 

.7 

.58 


TABLE I; Certifying Jacobian Non-Singularity. Maximum 
Value of 7 for different networks 


Fig. 4: 6 -bus network 


C. 14-bus system 
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Fig. 5: Cross-section of the true feasibility region and the 
inscribed certificate for IEEE 14-bus system. 


We repeat the same experiment with the IEEE 14 bus 
network. In this network, we find that a uniform bound on all 
injections, is too conservative. We thus use a heuristic to pick 
the relative scaling of the bounds on the injections at each bus. 
We create random voltage profiles satisfying the operational 
constraints, and compute the corresponding injection vectors. 
We then “fit” a box region to the resulting samples, and attempt 
to certify that this region (or a scaled version) is completely 
contained within the feasible set. 

We plot a projection of the feasible injection onto the pi—p 5 
space (the injections at buses 4 and 5) in figure The hgure 
shows that S indeed does touch the boundary of the feasible set 
and cannot be expanded further. The actual bounds established 
for each bus in the network are quire reasonable and are shown 
in table |II] 


D. Computation Time 

Table |I^ shows the size of the largest positive semidef- 
inite (PSD) constraint produced while solving the sparsity- 
exploiting moment relaxations of (0,(|^. This is usually the 
dominating factor in the computational effort of solving the 
moment relaxation. In | |27l , the authors show how a similar 
moment relaxation can be scaled upto a 300 bus network. In 
newer recent work, the authors have managed to scale this up 
to several thousand buses 128) . 

In terms of actual implementation, we use the convex 
optimization parser-solver CVX p9) p0| and MOSEK pi] 


as the underlying solver. Eor the 14-bus system, on a 2014 
Macbook Pro with a 2.6GHz Intel Core i7 processor, it takes 
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Bus 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

Pmin 

-0.2942 

-0.1455 

-0.3904 

-0.5806 

-O.I6I8 

-0.2446 

-0.1108 

-0.3683 

-0.1274 

-0.0939 

-0.0446 

-0.0953 

-0.0504 

Pmax 

0.5280 

0.1752 

0.7206 

0.5078 

0.3206 

0.3401 

0.1073 

0.3285 

0.3497 

0.1750 

0.1292 

0.2338 

0.1198 


Bus 

4 

5 

7 

9 

10 

11 

12 

13 

14 

9 min 

0.2802 

0.4454 

0.0253 

-0.0964 

0.0210 

0.0181 

-0.0086 

0.0020 

0.0066 

qmax 

0.7384 

0.8873 

0.1944 

0.1750 

0.2547 

0.1586 

0.1278 

0.2213 

0.1034 


TABLE II: Bounds on active and reactive power injections at each bus in the IEEE 14 bus network, guaranteeing existence of 
a feasible power flow solution (in per unit system with a 100 MVA base) 


20 seconds to solve and each instance of 0- Using the 
ideas from fig, we believe that the approach can be 

scaled to several thousand buses. 


Num of buses 

3 

6 

14 

30 

Size of PSD 

10 

55 

105 

190 


TABLE III: Certifying Jacobian Non-Singularity. Maximum 
Value of 7 for different networks 


VI. Discussion and Conclusions 


We have developed a novel computational approach for 
constructing geometrically simple regions of existence of 
feasible solution in the power injection space. Construction of 
these regions is based on semidefinite programming techniques 
which were recently proven to be both effective and scalable 


|, |27|. Our numerical experiments on small test networks 
have indicated that the constructed regions are tight, so any 
rescaling of a particular shape results in violation of constraints 
at least for some operating points. Regions that are certified 
this way have a very simple geometric shape of a polytope or 
an ellipse and can be naturally used in a number of security 
assessment or emergency control applications that require fast 
decision making. Two particularly suitable applications are the 
assessment of overload and voltage collapse risks in the pres¬ 
ence of renewable generation uncertainty and identification of 
optimal load shedding actions in emergency situations. 

There are several ways how the approach can be improved 
and extended. Eirst, the current version of the algorithm does 
not optimize with respect to the shape of the certified region. 
Identifying the largest ellipsoid that can be inscribed in a 
feasibility set may significantly improve the conservativeness 
of the approach. These techniques can be naturally extended 
to small-signal stability (analyzing Jacobians of the dynamical 
equations instead of the power flow equations). Integration of 
small-signal and transient stability for credible contingencies 
would complete the full characterization of the safe operation 
region and will become a powerful tool for system operators. 
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VII. Appendix 
A. Proof of theorem |///■ 7 | 

Let s S iS be arbitrarily chosen. Using assumption we 
know that 3s° S S and a valid voltage phasor vector U^such 
that s° = F (U°), H°P (U°; 7 ) > 0. We study the dynamical 
system 

with initial condition V (0) = V°. Then ^ (F (V) — s) = 
— {F {V) — s) so that 

P (0) — s = exp (—7) (F (U°) — s) => 

P (0) = 'S(l ~ exp (—7)) + exp (—7) G 5 (8) 

By ( |4al ), we know that Jp (U) is non-singular for every feasi¬ 
ble for every V such that F[°p (U; 7) > 0, {V) =0, so the 

dynamical system is well defined as long as H°p (U; 7) > 0 
(the constraint (U) = 0 is always maintained since 
the dynamics only changes , (|V|)pq). We start with 

jjop (V^'j > 0. Suppose that t* < 00 is the first time 
instant the dynamical system hits the operational boundary, 
ie, {V (t*)) = 0 for some i. The dynamical system is 
well-defined upto this point. By implication ( |4b| ), we know 
that since F (V (t*)) G S, H°p {V (t*)), (V {t*)) = 0, 
we must have H°p {V (t*)) > 0, which is a contradiction to 
F[°^ iy (t*)) = 0. Thus, the dynamical system always satisfies 


H°P (y y) > 0 , so that ([^ implies that it must converge to 
V such that F {V) = s, F[°p {V) > 0, so that s is strictly 
feasible. Since s G S was chosen arbitrarily, 5 is a domain of 
strict feasibility. 

B. Moment Relaxation 

In this section, we describe the moment relaxation of 
For any x G M”, let Poly* (x) denote the vector 
of all monomials of degree upto i \n x (with the first 
entry equal to the 0-degree monomial 1). For example, 
Poly^ (a;) = (l xi ... Xn x\ xix^ . Let 

V* = (z Re(U) Im(U))^. Let m be the size of 
Poly*' (V°). We define a moment vector y of the same size as 
Poly*' (U°) and the linear operator on the space of all degree-4 
polynomials in {V^). 

C m \ m 

^c,Poly4(n =Y,c,y, 

/ i-1 

We also define the localizing matrix for z = 1,2: Xi = 
Poly* (U'^) (Poly* (U*^))^. The moment relaxation of (|^ is 
given by the following feasibility problem; 

Find y such that (9a) 

Cy{JF{V)z{l U"^))=0 (9b) 

((l^i “ LT - -'^ 1 ) ^ 0 (9c) 

^y((|L|"-f-)Vi) =0,iGpv (9d) 

Cy{{\Vo-y)Xf)=Q,Cy{{z^z-l)Xi)=Q (9e) 
Cy{X2)FQ,yi = l (9f) 

Each constraint here is either a linear equality constraint or 
a linear matrix inequality in y. If ^ is feasible, then © is 
feasible as well; Let (z, V) be feasible for and define = 
(z Re(U) Im (U))^. Then, y = Poly^ (y*) is feasible for 
Q by construction. Thus, if (|^ is infeasible, then so is <§)■ 
The moment relaxation of problem Q can be formulated 
similarly. Suppose that 5 = {s : ||s —s'’|| < (ellipsoids 
and polytopes in s can be handled similarly). We define = 

(s Re(U) Im (U))^ and y, Poly*'(y^), Vi, similarly 

as before. We can then construct the moment relaxation of 0 ; 

Find y such that (10a) 

£y {yP (V) X,) = 0, Cy {{Fy {V) - S y) X ,) = 0, J = 1, . . . , k 

( 10 b) 

Cy{H;py)x,)=o,j = i,...,\H°p\,j^i (loc) 

£y{H:;yv)xf)h0,j = i,...,\Hy (lOd) 

£v (V 2 ) F 0, yi = 1, £y ((^<52 - ||s - ) Vi) F 0 

(lOe) 

Infeasibility of ( fTO) ! implies infeasibility of 0, using a similar 
argument as before. However, if ( fTO) ! is feasible, one can 
extract (U, s) from y. If this is feasible for 0 , then we know 
that S is not a DSF, since it contains a point that is not strictly 
feasible. 




